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The problem of corrections to Einstein's equations arising from averaging of inhomogeneities ("back- 
reaction") in the cosmological context, has gained considerable attention recently. We present re- 
sults of analysing cosmological perturbation theory in the framework of Zalaletdinov's fully covariant 
Macroscopic Gravity. We show that this framework can be adapted to the setting of cosmological 
perturbations in a manner which is free from gauge related ambiguities. We derive expressions for 
the backreaction which can be readily applied in any situation (not necessarily restricted to the 
linear perturbations considered here) where the metric can be brought to the perturbed FLRW 
form. In particular these expressions can be employed in toy models studying nonlinear structure 
formation, and possibly also in JV-body simulations. Additionally, we present results of example 
calculations which show that the backreaction remains negligible well into the matter dominated 
era. 



I. INTRODUCTION 

It is known [l| that in order to apply the equations of 
Einstein's General Relativity at the large length scales 
of interest in cosmology, it is necessary to first perform 
a smoothing or averaging operation, which will gener- 
ate nontrivial corrections in these nonlinear equations. 
There is a debate in the literature concerning two basic 
questions regarding this smoothing operation : (a) How 
does one obtain consistent and observationally relevant 
variables associated with the corrections due to averag- 
ing or "backreaction" , and (b) what is the magnitude of 
this backreaction? 

This discussion has had a long history 0, H, 0], al- 
though in recent times attention has been focussed on 
two promising candidates for a consistent nonperturba- 
tive averaging procedure, jiamely the spatial averaging 
of scalars due to Buchert [j| and the covariant approach 
due to Zalaletdinov [f| 0] ■ Perhaps owing to its appeal- 
ing simplicity of implementation, the approach due to 
Buchert has attracted a significant amount of attention 
both in the context of cosmological perturbation theory 
a [ 9. flfj| and with regards fully nonlinear calculations 
111 Ha . Il3| , although very interesting results have been 
derived using Zalaletdinov's covariant procedure as well 
0, EH • The most fascinating aspect of these studies has 
been the possibility that the phenomenon of Dark Energy 
[1, [lU and also Dark Matter [l5| might be attributed to 
the backreaction from averaging (see Ref. [l6| for a re- 
cent review of the subject). 

The physical relevance of the averaged variables de- 
fined in Buchert 's formalism, has been questioned in the 
literature [ItJ ■ It has been argued [l?], [3 that effects 
of inhomogeneities which only perturbatively affect the 
(on average homogeneous) metric of the universe, cannot 
lead to effects large enough to account for the inferred 



acceleration of the universe from, say supernovae obser- 
vations [19| . This has been countered by the argument 
[HI that perturbation theory may break down during the 
epoch of fully nonlinear structure formation. Recently it 
was shown [20( in the context of a specific model of spher- 
ical collapse, that an explicit coordinate transformation 
could be found, which brought the metric to the per- 
turbed Friedmann-Lemaitre-Robertson- Walker (FLRW) 
form, 
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-(1 + 2<p)dr 2 + a 2 (r)(l - 2iJ;)dx 2 
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satisfying all conditions required for a perturbation for- 
malism in if and ip to hold, even in the regime of fully 
nonlinear collapse. 

In such a context, it is important to have a consis- 
tent formalism, free from issues such as gauge artifacts, 
in which one has derived expressions for the backreac- 
tion which can be applied in a straightforward manner 
to any model in which the metric can be brought to the 
form (flj. Since the Buchert framework by construction 
is best adapted to coordinates comoving with the mat- 
ter, most applications using perturbation theory in this 
framework have focussed on the synchronous and comov- 
ing gauge, although recently Behrend et al. [lfj have 
also performed calculations in the conformal Newtonian 
gauge. However, as pointed out elsewhere [2(|, calcula- 
tions using the Buchert framework in perturbation theory 
necessarily face the ambiguity of dealing with two scale 
factors : one is the scale factor with which the perturbed 
FLRW metric is defined, and the other is the volume av- 
eraged scale factor defined by Buchert [f|, and it is not 
clear which of these scale factors is the observationally 
relevant one. 

In this paper we shall deal with the covariant frame- 
work of Zalaletdinov's Macroscopic Gravity (MG), and 
its spatial averaging limit proposed in [2l|. We will see 
that here, one has a well defined averaged metric with 
a corresponding uniquely defined scale factor. Further, 
we will argue that the structure of MG allows us to com- 
pletely specify a consistent, observationally relevant aver- 



2 



aging operation adapted for perturbations to the FLRW 
geometry. The final expressions obtained for the backre- 
action can be applied directly (or with some straightfor- 
ward modifications), to any situation in which the metric 
can be brought to the form ([I]) or a suitable generalisa- 
tion thereof. 

The organisation of this paper is as follows : Sec. |TT] 
presents a brief recap of the basic ideas underlying MG 
and some useful expressions from standard cosmological 
perturbation theory (PT), together with a proposal for 
practically estimating the effect of the backreaction on 
the time evolution of the FLRW scale factor. Sec. IIIII 
presents details of the MG averaging procedure adapted 
to cosmological PT, including general expressions for the 
leading order backreaction terms, with a discussion of 
gauge related issues and the definition of the averaging 
operator. The heart of the paper is in Sec. IIV1 where 
we derive final expressions for the backreaction, both 
in real space and Fourier space, which can be directly 
utilised in model calculations. While these expressions 
use a few simplifying restrictions, these can be lifted if 
necessary in a completely straightforward manner. Sec. 
IVl contains example calculations in first order PT, which 
show that the magnitude of the backreaction is, as ex- 
pected, negligible compared to the homogeneous energy 
density of matter in the radiation dominated era and for 
a significant part of the matter dominated era. We con- 
clude in Sec. I VII with some final comments. Throughout 
the paper, lower case Latin indices a,b,c... will refer to 
spacetime indices 0, 1, 2, 3, and upper case Latin indices 
A,B,C,... to spatial indices 1,2,3. The speed of light 
c is set to unity, and a prime refers to a derivative with 
respect to conformal time unless stated otherwise. 



calculations is 



W a b (x',x) 



dx a 



<9Xyp 



dx b 
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where Xy P refers to a coordinate system in which the 
metric determinant is constant. The class of such co- 
ordinate systems is called volume preserving (VP). Not 
surprisingly, the entire formalism of MG simplifies con- 
siderably in VP coordinate systems, in which the coordi- 
nation bivector reduces to the Kronecker delta 



m (x', x ) 



VPC 



(3) 



We will see later that VPC systems play an important 
role in consistently setting up cosmological perturbation 
theory in the context of MG. 

The average of a tensor P% over a finite spacetime do- 
main S is given by 



P a h (x) = ( Pi )(x) 



(4) 



where Pu(x', x) is the bilocal extension of P% defined as 



P a b (x',x) = W?,(x,x')P),(x')Wl {x',x) 



(5) 



The averaging operation when appropriately applied to 
the connection T bc on Ai, gives an averaged connection 
rg c which is taken to be the connection on an averaged 
manifold M. In other words, the connection on M. sat- 
isfies the condition 



II. COVARIANT MACROSCOPIC GRAVITY 
(MG) AND COSMOLOGY 

A. Recap of MG formalism 



In this section we present a rapid overview of the 
generally covariant averaging formalism of Macroscopic 
Gravity (MG), developed by Zalaletdinov and coworkers 
@> Bi J22I [23( 1 . and its spatial averaging limit proposed 
in [21] . The reader is referred to these papers for more 
details. 

The MG formalism uses a bilocal operator W* (x',x), 
called the coordination bivector, to define a covariant 
spacetime averaging operation on tensors in a spacetime 
manifold A4. [The prime here is being used to distinguish 
two different spacetime points, and must not be confused 
with a derivative with respect to conformal time. ] The 
various properties that this operator must satisfyand a 
proof of its existence can be found in Refs. @, @, 
The form of the coordination bivector used in all MG 
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The metric G a b associated with the averaged connection 
can be assumed to be the average of the inhomogeneous 
metric g ab on M, i.e. G a b — 9ab- (We will see later that 
in the perturbative setting this amounts to a very natural 
condition on the perturbations.) Averaging the Einstein 
equations on M leads to the equations satisfied by the 
averaged metric, which can be written as 



El 



8nG N Tu + i3rav) T a _ 



(7) 



Here E% is the Einstein tensor constructed from the met- 
ric Gab and its inverse G ab , T b is the averaged energy- 
momentum tensor, Gat is Newton's gravitational con- 
stant, and (9™v) T a 

is a tensorial correlation object which 
acts like an effective gravitational energy-momentum ten- 
sor. For brevity we shall omit its detailed definition, re- 
ferring the reader to Ref. 0, 0] , and only note that its 
broad structure can be symbolically represented as 



(grav)rp 



(8) 
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where T symbolically denotes the bilocal extension of the 
Christoffel connection on Ai. In general the total energy- 
momentum tensor (8ttGnT^ + (grav)rpa^ j g CO variantly 
conserved, 



^G N T a b +^rav) T ay a =Q 



(9) 



with the semicolon denoting the covariant derivative with 
respect to the averaged geometry (f bc and G a b)- 

In Ref. [2l[ it was argued that in the cosmological con- 
text, it is essential to consider a spatial averaging limit 
of the covariant averaging used in MG. The simplest way 
to see this is to note that the homogeneous and isotropic 
FLRW spacetime must be left invariant under the aver- 
aging operation, and this is only possible if the averaging 
is tuned to the uniquely defined spatial slices of constant 
curvature in the FLRW spacetime. A related statement 
also applies to perturbation theory in general. The appli- 
cation of the MG formalism to perturbation theory has 
been discussed in Refs. 0, HH, which also compare the 
results of MG in the perturbative scenario, with those 
obtained by Isaacson [2J| in the study of the short wave- 
length limit of gravitational waves. A crucial observation 
[23j is that what one normally refers to as the "back- 
ground" in perturbation theory, must be defined by an 
averaging procedure. In particular, the background must 
remain invariant under the averaging operation. In Ref. 
[2l| no a priori assumption was made about the form 
of the inhomogeneous metric, and hence certain assump- 
tions had to be made about the choice of spatial slicing 
used to define the spatial averaging. In the case at hand, 
namely when the inhomogeneous metric is taken to be a 
perturbation around the FLRW metric, we will see that 
the situation simplifies to some extent and a consistent 
spatial averaging operation can be identified. 

Throughout this paper we will assume that the met- 
ric of the universe is a perturbation around the FLRW 
metric given by 



ds 2 = a 2 (r/) (-dr/ 2 +j AB dx A dx B ) 



(10) 



Here a is the scale factor and r\ is the conformal time 
coordinate related to cosmic time r by the differential 
relation 



dr = a(rj)dr/ 



(11) 



In Eqn. (JTDJ) we have allowed the spatial metric to have 
the general form a 2 jab where 7,45 is the metric of a 
3-space of constant curvature. For the calculations in 
this paper we will assume a flat FLRW background in 
coordinates such that 7,43 = Sab, however for future 
reference we shall present certain expressions in terms of 
the more general spatial metric. 

Assuming the averaged metric to have the form ()10D . 
and the averaged energy-momentum tensor to have the 
form 



b 



{ P + p)v a v b +p5 a b 



b ' 



(12) 



where v a is the timelike 4-vector which defines the ho- 
mogeneous spatial slices of the FLRW spacetime, and 
p{r) and p(r) are respectively the homogeneous energy 
density and pressure corresponding to the averaged mat- 
ter distribution, the modified cosmological equations ob- 
tained by averaging a general (i.e. not necessarily per- 
turbed FLRW) spacetime, can be shown to reduce to the 
following (see Eqns. 87 of Ref. [U) 



1 da\ _ 8ttG n 1 
a dr J 3^6 



p(D +lS (D 



(13a) 



1 d 2 a AttG n 1 
ad^ = -— {P + 3p)+ 3 



■p(l) + 7?(2) +< §(2) 



(13b) 



where the combinations (P (1) + <S (1) ) and ("P (1) +V ( - 2) + 
S^) are generally covariant scalars defined by the rela- 
tions (see Eqns. 88 of Ref. [2l[ , with f = a) 
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(14b) 
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5 (2) = i< r^f% ) + (g JK )( r° JQ r£ A ) , (i4d) 



where we have defined 



n 



1 da 
a dr] 



a 
a 



(15) 



and accounted for the fact that in general, the average 
of the inverse inhomogeneous metric, need not equal the 
inverse metric G ab . However, we will soon see that in the 
perturbative setting we can in fact set these two tensors 
to be equal via a very natural condition on the perturba- 
tions. The index in Eqns. pT|) refers to the conformal 
time rj. The averaging in Eqns. (|14p is assumed to be a 
spatial averaging in an unspecified spatial slicing in the 
inhomogeneous manifold Ai; in Sec. IIIII we will specify 
the averaging procedure more exactly. 

In addition, the following "cross-correlation" con- 
straints must also be satisfied by the inhomogeneities (see 
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Eqns. 89 of Ref. [2J) 

+ {g 



~JK \ I r r B 
\ 1 ,/S 1 AK 



= (16a) 



~JK 



A r B 
JB L OK 



) - ( r%r» K ) 



= (16b) 



/ -pm \ / -pA pm 

\ 1 BO 1 Om / \ 1 ml) 1 OS 



1 Jm 1 A'S 



) ( ^JB^^Km ) 



°B 



3 V 



4ft 2 



(16c) 



where the lower case index m in the last equation runs 
over all spacetime indices 0,1,2,3. 



B. Cosmological Perturbations and Gauge 
Transformations 



ip,ip,u> and x each correspond to one scalar degree of free- 
dom, the transverse 3- vectors Hi a and \A each correspond 
to two functional degrees of freedom, and the transverse 
tracefree 3-tensor xab corresponds also to two functional 
degrees of freedom. This totals to 10 degrees of freedom, 
of which 4 are coordinate degrees of freedom which can 
be arbitrarily fixed, which is what one means by a gauge 
choice. For example, the conformal Newtonian or lon- 
gitudinal or Poisson gauge [2a . l27l |28| is defined by the 
conditions 

uj = Q = X ; x A = o. (21) 

For the metric (fl~7|) we have at first order, 

yj-dct g = a 4 (r/) (1 + tp - Sip) . (22) 

The inverse of metric (fl7j). correct to first order, has the 
components 
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(l-2p) ; <?' 



OA 



1 
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AB 



a a 
-1 ((1 + 2^)7^ - X AB ) 



(23) 



Denoting 7i = (a'/a), the prime denoting a derivative 
with respect to conformal time r/, the first order accurate 
Christoffel symbols are 



For ready reference, in this subsection we present ex- 
pressions for the metric, its inverse, and the Christoffel 
connection in first order cosmological PT, in an arbi- 
trary, unfixed gau ge. The notation we use is similar to 
that used in Ref. [2a]. We will also give expressions for 
the first order gauge transformations of the perturbation 
functions (see e.g. Ref. (26|). 

The first order perturbed FLRW metric in an arbitrary 
gauge can be written as 



1 oo 



H + <p' ; r° 0A = d Al p + HLO A 



r A = d A p + uj Ai + Hu A , 



r° AB = (n - v/ - m{p + V)) iab - y ( A^ B) 

+ \x'AB + KXAB , 



ds 2 = a 2 (77) [-(1 + 2ip)dn 2 + 2uo A dx A dn+ 

((1 - 2iP) lAB + xab) dx A dx B ] . (17) 



OS 



(H - if) S 



v Al 

:Xb ■ 



The functions tp and tp are scalars under spatial coordi- 
nate transformations. The functions uj a and xab can be 
decomposed as follows 



UA = 9a^> + &A 



XAB = D AB X + 2V (j4 Xs) + XAB , 

(18) 



where the parentheses indicate symmetrization; D A b is 
the tracefree second derivative defined by 



D AB = V A V B - (l/3) 7A sV 2 



with V a the covariant spatial derivative compatible with 
7as; and oj A , xa and xab satisfy 

V a lo a = Q = VaX A ; V A Xs=0 = X^, (20) 

where spatial indices are raised and lowered using j A b 
and its inverse j AB . From their definitions it is clear that 



BC 



(3)f A 
1 BC 



(S A d c tP + 6^d B ^ - lB C d A ^) 



Hu A lBC + - (VcXb 



VsXc " V^xsc) 



(24) 



where ^Tgc denotes the Christoffel connection associ- 
ated with the homogeneous 3-metric 7,4s- 

Gauge transformations : While the concept of 
gauge transformations can be described in a rather so- 
phisticated language using pullback operators between 
manifolds 26], for our purposes it suffices to implement 
a gauge transformation using the simpler notion of an in- 
finitesimal coordinate transformation (also known as the 
"passive" point of view) [2{|. Hence, denoting the co- 
ordinates and perturbation functions in the new gauge 
with a tilde (i.e. x a , tp, lu a , and so on), we have 



x a +£ a {x) ; x a 



(25) 
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where the infinitesimal 4- vector £ a can be decomposed as 

t a ={e,Z A ) = (a,d A l3 + d A ) , (26) 

where a and /3 are scalars and d A is a transverse 3-vector 
satisfying Vag^ 4 = 0. 

It is then easy to show that if this transformation is 
assumed to change the metric (TT71) by changing only the 
perturbation functions but leaving the background intact 
(a so-called "steady" coordinate transformation), then 
the old perturbations and the new are related by [26j 

ip = (p + a' + Ha , 

ui = u) — a + f3' , 
Cj a = h A + d Al , 
X = X + 2/3 , 

x A = x A + d A , 

XAB = XAB ■ (27) 



The last equality shows that the transverse tracefree ten- 
sor perturbations are gauge invariant. They correspond 
to gravitational waves. 



Time evolution of the background : An 
iterative approach 



artificially set C to zero and obtain and ip^ using 
the standard approach (see e.g. Ref. [20]). Clearly, since 
the "true" background (say a*) satisfies Eqns. (H2J) with 
a nonzero C, we have in general =/= a*. Now, using 
the solution (p^ ', we can calculate the zeroth iteration 
correlation objects by applying the prescription to 
be developed later in this paper. As a first correction to 
the solution , we now solve for a new background a*- 1 -* , 
with the known functions acting as sources in Eqns. 
(TL3")) . This first iteration will then yield a solution (p^ for 
the inhomogencitics, and hence a new set of correlation 
terms and this procedure can be repeatedly applied. 
[See however the first paper in Ref. L4| for an alternative 
approach exclusively using averaged quantities in solving 
the full MG equations.] Pictorially, 



(o) 



(i) 



(28) 



As for convergence, if perturbation theory is in fact a 
good approximation to the real universe, then one can 
expect that the correlation terms will tend to be small 
compared to other background objects, and will therefore 
not affect the background significantly at each iteration, 
leading to rapid convergence. On the other hand, if the 
correlation terms are large, this procedure may not con- 
verge and one might expect a breakdown of the pertur- 
bative picture itself 31]. We will see that in the linear 
regime of cosmological perturbation theory, the correla- 
tion terms do in fact remain negligibly small. 



Before we move on to deriving formulae for the corre- 
lation terms (|14[) in terms of perturbation functions in 
the metric, there is one issue which merits discussion. 
The cosmological perturbation setting, together with the 
paradigm of averaging, presents us with a rather pecu- 
liar situation. On the one hand, the time evolution of the 
scale factor is needed in order to solve the equations satis- 
fied by the perturbations. Indeed, the standard practice 
is to fix the time evolution of the background once and 
for all, and to use this in solving for the evolution of the 
perturbations. On the other hand, the evolution of the 
perturbations (i.e. - the inhomogeneities) is needed to 
compute the correlation terms appearing in Eqns. ([13]) . 
Until these terms are known, the behaviour with time of 
the scale factor cannot be determined; and until we know 
the scale factor as a function of time, we cannot solve for 
the perturbations. Note that this is a generic feature 
independent of all details of the averaging procedure. 

It would appear therefore, that we have reached an 
impasse. To clear this hurdle, one can try the following 
iterative approach : Symbolically denote the background 
as a, the inhomogeneities as ip, and the correlation ob- 
jects as C . Note that a, ip and C all refer to functions 
of time. We start with a chosen background, say a stan- 
dard flat FLRW background with radiation, baryons and 
cold dark matter (CDM), and solve for the perturbations 
in the usual way, without accounting for the correlation 
terms C. In other words, for this "zeroth iteration" , we 



III. THE AVERAGING OPERATION AND 
GAUGE RELATED ISSUES 

In this section, we will describe the details of the MG 
(spatial) averaging procedure adapted to the setting of 
cosmological PT. 

A. Volume Preserving (VP) Gauges and the 
Correlation Scalars 

It will greatly simplify the discussion if we start with 
symbolic calculations which allow us to see the broad 
structure of the objects we are after. Since the corre- 
lation objects in Eqns. (fl~3|) depend only on derivatives 
of the metric, we will primarily deal with metric fluc- 
tuations; matter perturbations will only come into play 
when solving for the actual dynamics of the system. Be- 
fore dealing with the issue of which gauge to choose in 
order to set the condition ([B]) , we will show that irrespec- 
tive of this choice, the leading order contribution to the 
correlations requires knowledge of only first order per- 
turbation functions. 

We will use the following symbolic notation : 

• Inhomogeneous connection: F 

• FLRW connection: Tp 
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• Perturbation in the connection : <5r = T — Tp = 

sri 1 ) + srw + ... 

• Coordination bivector : W = 1+SW = 1 + W (1) + 
6W& + ... 

• Bilocal extension of the connection : T 

• Inhomogencous part of the bilocal extension of the 
connection : ST = f-T F = ST W + of (2) + ... 

• Correlation object : (s™ u )T 

The integer superscripts denote the order of perturba- 
tion. The form of the coordination bivector arises from 
the fact that in perturbation theory, in the spatial av- 
eraging limit, a transformation from an arbitrary gauge 
to a VP one can be achieved by an infinitesimal coordi- 
nate transformation. By a VP gauge we mean a gauge 
in which the metric determinant is independent of the 
spatial coordinates to the relevant order in PT, but may 
be a function of time. It can be shown that such a func- 
tion of time (which will typically be some power of the 
scale factor), is completely consistent with all definitions 
and requirements of MG in the spatial averaging limit. 
An easy way of seeing this is to note that in any av- 
eraged quantity, the metric determinant appears in two 
integrals, one in the numerator and the other in the de- 
nominator (which gives the normalising volume). In the 
"thin time slicing" approximation we are using to define 
the averaging, any overall time dependent factor in the 
metric determinant therefore cancels out. Also, a fully 
volume preserving coordinate system can clearly be ob- 
tained from any VP gauge as defined above, by a suitable 
rescaling of the time coordinate. It is not hard to show 
that in the thin time slicing approximation, this gives the 
same coordination bivector (x',x) as the VP gauge 
definition above. 

To see that first order perturbations are sufficient to 
calculate (srav)j^ ^ Q i ea( j m g order, we only have to note 
that the background connection satisfies 

(T F ) = T F , (29) 

and that the structure of O^T is given by Eqn. ©. 

(grav) T then reduces to 

(grav) T ^^^j^ff^ ( 30 ) 

which is exact. Clearly, the correlation is quadratic in the 

perturbation as expected, and hence to leading order, i5r 

~(l) 

above can be replaced by ST 

Eqns. (|2"9"|) and (|3T))) treat the averaging operation at a 
conceptual level only. To make progress however, we also 
need to prescribe how to practically impose the averaging 
assumption 



in any given perturbative context. This requires some 
discussion since, for example, the bilocal extension of the 
connection T has the structure 



T = TW + W- 1 (d + d')W , 



(32) 



where d is a derivative at x and & a derivative at x' . 
[The reader is referred to Ref. [1, 0] for the detailed ex- 
pression. Suffice it to note that this structure ensures 
that the averaged connection has the correct transfor- 
mation properties.] The actual MG averaging operation 
in general is therefore a rather involved procedure. Ad- 
ditionally, it is also necessary to address certain gauge 
related issues. 

To clarify the situation, let us start with a fictitious 
setting in which the geometry has exactly the flat FLRW 
form, with no physical perturbations. Clearly, if we work 
in the standard comoving coordinates in which the metric 
jab of Eqn. (fT0|) is simply jab — Sab, then since these 
coordinates are volume preserving in the sense described 
above, the coordination bivector becomes trivial. The 
averaging involves a simple integration over 3-space, and 
we can easily see that Eqn. (|29[) is explicitly recovered. 

Now suppose that we perform an infinitesimal coordi- 
nate transformation, after imposing Eqn. (129p . Since the 
averaging operation is covariant, then from the point of 
view of a general coordinate transformation, both sides 
of Eqn. (f2"9"| will be affected in the same way. However, 
suppose that we had performed the transformation before 
imposing Eqn. (|29|) . In the language of cosmological PT, 
we would then be dealing with some "pure gauge" per- 
turbations around the fixed, spatially homogeneous back- 
ground. If we did not know that these perturbations were 
pure gauge, we might naively construct the nontrivial co- 
ordination bivector for this metric, compute the bilocal 
extension of the connection according to Eqn. (f3"2"| and 
try to impose Eqn. pip . This would be incorrect since 
these perturbations were arbitrarily generated and need 
not average to zero (for example they could be positive 
definite functions). In order to maintain consistency, it 
is then necessary to ensure in practice that the averaging 
condition (|3ip is applied only to gauge invariant fluctua- 
tions (which is rather obvious in hindsight). 

There is another problem associated with the structure 
of the coordination bivector, even when there are real, 
gauge invariant inhomogeneities present. Note from Eqn. 
© that the coordination bivector has the structure 



W = 



dx 



dx\ 



dx\ 



Ox 



(33) 



where x denotes the coordinates we are working in and 
xy a set of VPCs. In perturbation theory (in the spatial 
averaging limit) we will have, at leading order, 



x = x v - £ ; x v = x + £ , 



(34) 



where £ symbolically denotes an infinitesimal 4-vector 
defining the transformation, and hence 



r } = r 



f i.e. 



(ST)=0, 



(31) 



(dx v )/(dx) = 1 + d£ 



(35) 
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and so on, which gives us 

W = 1 - {dC)\x> + (dOU + . . . = 1 + SW™ + ... (36) 

Now when we compute a quantity such as ( TpSW^ ) 
which appears in the expression (|32p for ( T } , we will 
be left with a fluctuating (x-dependent) term of the form 
IV(( <9£ )—d£), where x denotes the 3 spatial coordinates. 
Hence if we try to impose Eqn. (|3ip we will ultimately 
be left with equations of the type 

( / ){x) - f{x) = (37) 

for some functions derived from the inhomogeneities 
which we have collectively denoted /. In other words, 
consistency would seem to demand that the inhomo- 
geneities vanish in this coordinate system, which is of 
course not desirable. 

It therefore appears that we are forced to impose Eqn. 
pip in a volume preserving gauge, since by definition, 
only in such a gauge will we have W — 1 exactly. We 
emphasize that this is a purely practical aspect related 
to defining the averaging operation, and is completely 
decoupled from, e.g. the choice of gauge made when 
studying the time evolution of perturbations. We are 
in no way breaking the usual notion of gauge invariance 
by choosing an averaging operator. The conditions Eqn. 
PT|) now reduce to the form 

( fvpc }(x)=0 (38) 

which are far more natural than Eqn. (|37p . The aver- 
aging condition is now unambiguous, but depends on a 
choice of the VP gauge which defines the averaging op- 
eration, an issue we shall discuss in the next subsection. 
For now, all we can assert is that this VP gauge must 
be such that in the absence of gauge invariant fluctua- 
tions, it must reduce to the standard comoving (volume 
preserving) coordinates of the background geometry as 
in Eqn. (jTUjl . This of course is simply the statement that 
the VP gauge must be well defined and must not contain 
any residual degrees of freedom. 

The averaging operation now takes on an almost trivial 
form - to leading order it is easy to show that for any 
quantity f(rj, x) (with or without indices), the average of 
/ in a VP gauge in the spatial averaging limit, is given 
by 

(f)( V ,S) = ^[ d 3 yf(n,y), (39) 

V L Jv(x) 

where the integral is over a spatial domain V(x) with 
a constant volume V^. The spatial coordinates are the 
comoving coordinates of the background metric, and at 
leading order the boundaries of V(x) can be specified in 
a straightforward manner as, e.g., 

V(x) = {y | x A - L/2 <y A <x A + L/2, A =1,2, 3.} , 

(40) 



where L is a comoving scale over which the averaging 
is performed (in which case Vl — L 3 ). The averaging 
definition can be written more compactly in terms of a 
window function Wj, (x, y) as 

( / )[r},x) = J d 3 yW L {x,y)f{n,y), 

J d 3 yW L {x,y) = l, (41) 

where Wl {x, y) vanishes everywhere except in the region 
V{x), with the integrals now being over all space. This 
expression will come in handy when working in Fourier 
space, as we shall do in later sections. 

A couple of comments are in order at this stage. 
Firstly, we have not specified the magnitude of the aver- 
aging scale L. The general philosophy is that this scale 
must be large enough that a single averaging domain en- 
compasses several realisations of the random inhomoge- 
neous fluctuations, and small enough that the observable 
universe contains a large number of averaging domains. 
However, as we will show later in Sec. IIV1 if one is ulti- 
mately interested in quantities which are formally aver- 
aged over an ensemble of realisations of the universe (as 
is usually done in interpreting observations), then the 
actual value of the averaging scale becomes irrelevant. 

This brings us to the second issue. The above dis- 
cussion is valid only in the situation where there are no 
fluctuations at arbitrarily large length scales, since in the 
presence of such fluctuations the averaging condition (f3"T) 
loses meaning (in such a situation it would be impossible 
to isolate the background from the perturbation by an 
averaging operation on any finite length scale). Indeed, 
we shall see a manifestation of this restriction in Sec. IIV1 
where the correlation scalars will be seen to diverge in 
the presence of a nonzero amplitude at arbitrarily large 
scales, of the power spectrum of metric fluctuations. 

We will end this subsection by explicitly writing out 
the averaging condition in an "unfixed VP" gauge, to 
be defined below, and also writing the correlation terms 
appearing in Eqn. (fT"3"|) . in this gauge. As we can see 
from Eqn. (f2"2")h the basic condition to be satisfied by a 
VP gauge is 

ip = . (42) 

Hereafter, all VP gauge quantities will be denoted using 
a tilde. This should not be confused with the similar no- 
tation that was used for the bilocal extension in Sec. Ill A[ 
which will not be needed in the rest of the paper, (p and 
ip are the scalar potentials appearing in the perturbed 
FLRW metric (JT7J). The single condition (O leaves 3 
degrees of freedom to be fixed, in order to completely 
specify the VP gauge one is working with. The MG for- 
malism by itself does not prescribe a method to choose 
a particular VPC system; in fact this freedom of choice 
of VPCs is an inherent part of the formalism. We shall 
return to this issue in the next subsection. For now we 
define the "unfixed VP (uVP) gauge" by the single re- 
quirement (|42p . with 3 unfixed degrees of freedom, and 
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present the expressions for the averaging condition and 
the correlation scalars, with this choice. 

It is straightforward to determine the consequences of 
requiring Eqn. ^ to hold, with the right hand side cor- 
responding to the FLRW connection in conformal coor- 
dinates, and remembering that the coordination bivector 
(in the spatial averaging limit) is now trivial (see Eqn. 
©). Together with some additional reasonable require- 
ments, namely 



V 2 s ) = = ( V 2 d A s 



(43) 



for any scalar s(r), x), the averaging condition in the uVP 
gauge reduces to 

{u> A ) = = {u>' A ) ; (x'ab) = 0, 

( VcXb ) + ( V B Xc ) ~ ( V a xbc ) = 0, 

( V^ B ) - ( V B u> A }=n(xA B ), (44) 

where we have used the expressions in Eqn. (|24[l with 
the uVP condition f4"2")) . We will also make the additional 
reasonable requirement that 



( XAB ) =0, 



(45) 



using which it is easy to see that the perturbed FLRW 
metric (IT7|) and its inverse (|2"3")) , in the uVP gauge, both 
on averaging reduce to their respective homogeneous 
counterparts, namely 



\ (FLRW) 

9ab )=g ab 



I n ab \ _ ab 

\ y I — "(FLRW) 



(46) 



Using these results, the expressions (|T4")) simplify to give, 
in the uVP gauge, 



l r 

~n2 



-\(XabX AB ') 



(47a) 



S (i) 



-10( d A ^d A ^ ) - 2( d A ^BX AB ) 



V B X AC (2V aX bc ~ VbXac) ) 



(47b) 



1 



V (l) +V (2) = J_ 



6( (^') 2 ) - 24H( ) 



-(iP'V 2 u>) + ~(x'abV A u) B ) 



-( Xab (x 



-ABi 



2H X 



~AB\ 



(47c) 



S {2) = -2 
a z 



3( Co A 'd A i> ) + H( Cj a u>' a ) 



(47d) 



where square brackets denote antisymmetrization. 



B. Choice of VP Gauge 

In this subsection we will prescribe a choice for the VP 
gauge which defines the averaging operation. In general, 
the class of volume preserving coordinate systems for any 
spacetime, is very large (see Ref. [22] for a detailed char- 
acterisation). We have so far managed to pare it down 
by requiring that the VP gauge we choose should re- 
duce to the standard FLRW coordinates in the absence 
of fluctuations. It turns out to be somewhat difficult to 
go beyond this step, since there does not appear to be 
any unambiguously clear guiding principle governing this 
choice. We will therefore motivate a choice for the VP 
gauge based on certain details of cosmological PT which 
one knows from the standard treatments of the subject. 

In particular, we shall make use of certain nice proper- 
ties of the conformal Newtonian or longitudinal or Pois- 
son gauge, which is defined by the conditions (|2"Tj) [261.I27T] 
(henceforth we shall refer to this gauge as the cN gauge 
for short). Since this gauge is well defined and has no 
residual degrees of freedom, all the nonzero perturbation 
functions in the cN gauge, namely tp, ip, Cj a and xab 
in the notation of Sec. Ill Bl are equal to gauge invari- 
ant objects. This is trivially true for xab, as seen in 
the last equation in (|27p . For the rest, note that in any 
arbitrary unfixed gauge, the following combinations are 
gauge invariant at first order 



$s =ip+-d„ 
a 



a I w- -X 



1 ./ 



Va = &a - Xa , 



(48) 



which can be easily checked using Eqns. (|2"T|) , and in the 
cN gauge, u>, \ and xa all vanish. Here and ^>b are 
the Bardeen potentials [32j (upto a sign), and \&s in par- 
ticular has the physical interpretation of giving the gauge 
invariant curvature perturbation, which is the quantity on 
which initial conditions are imposed post inflation [33j . 

Additionally, it is also known that the cN gauge for the 
metric remains stable even during structure formation, 
when matter inhomogeneities have become completely 
nonlinear. (See Ref. [l7J for an intuitive description of 
why this is so, and Ref. [2p| for an explicit demonstra- 
tion in a toy model of structure formation.) We believe 
that this is a strong argument in favour of using the cN 
gauge to define a VP gauge which will then define the av- 
eraging operation in the perturbative context. This will 
ensure that this "truncated" averaging operation, defined 
for first order PT, will remain valid at leading order even 
during the nonlinear epochs of structure formation. 

To implement this in practice, consider a transforma- 
tion from the cN gauge to the uVP gauge defined by Eqn. 
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(|42j) . The transformation equations (p7|) reduce to 
a' + ma + V 2 /3 = tp - 3tp , 
■iJj — — p — a' — ?fo , 

uj = a — (3' , 
^ = Cj a - d Al , 
X = -2/3, 

XAB = XAB 



(49) 



Recall that to completely specify a VP gauge, we need to 
fix 3 degrees of freedom in the uVP gauge. Our require- 
ment regarding the "well defined" -ness of the VP gauge, 
forces us to set d A — 0, and to choose a and (3 such that 
they vanish in the case where ip = = ip. 

This has fixed 2 degrees of freedom, in addition to the 
condition (j42f which is just the definition of the uVP 
gauge, and has hence not yielded a uniquely specified VP 
gauge. To do this, we shall make the following additional 
requirement. Since we arc dealing with a spatial averag- 
ing, it seems reasonable to require that the VP gauge 
being used to define the averaging, should be "as close as 
possible" to the cN gauge in terms of time slicing, and 
for this reason we shall set the function a to zero. 

To summarize, the VP gauge chosen is defined in terms 
of the gauge transformation functions £ a = (a, d A fi + d A ) 
between the cN gauge and the VP gauge, by the following 
relations 



and 



a = = d A , 

ip = 3ip = p , 

V 2 /3 = ip - Zip , 

Q = -p' ; x=~2/3, 

Xa = 0, 

Qj a =lo a ; Xab = XAB , 



(50) 

(51a) 
(51b) 
(51c) 

(51d) 
(51e) 



where the function is restricted not to contain any 
nontrivial solution of the homogeneous (Laplace) equa- 
tion V 2 /3 = 0. 

Having made this choice for the VP gauge, we are now 
assured that all averaged quantities which we compute 
are gauge invariant : our choice ensures that the averag- 
ing procedure does not introduce any pure gauge modes, 
and the philosophy of "steady" coordinate transforma- 
tions ensures that all background objects are, by assump- 
tion, unaffected by gauge transformations. In particular, 
the correlation objects in Eqns. (fT4")l arc all gauge in- 
variant. This is different from the gauge invariance con- 
ditions derived in the first paper of Ref. [1], where the 
background was also taken to change under gauge trans- 
formations at second order in the perturbations. It is at 
present not clear how these results are related to ours. 



Note that all these arguments are valid at first order 
in PT, which is sufficient for our present purposes. A 
consistent treatment at second order would require more 
work, although as long as one is interested only in the 
leading order effect, these arguments are expected to go 
through. 



IV. THE CORRELATION SCALARS 

With the VP gauge choice defined by Eqns. (fijTj) , it 
is straightforward to rewrite the correlation objects in 
Eqns. (|47f (which are in the uVP gauge) in terms of 
the perturbation functions in the cN gauge. We will re- 
strict the subsequent calculations in this paper to the 
case where there are no transverse vector perturbations, 
i.e., 



&A — : 



(52) 



in the cN gauge. This is a reasonable choice since such 
vector perturbations, even if they are excited in the initial 
conditions, decay rapidly and do not source the other 
perturbations at first order |30|]. 

In addition, for simplicity (and to keep this paper con- 
cise) we will choose to ignore the gauge invariant tensor 
perturbations as well, 



Xab = . 



(53) 



It will be an interesting excercise to account for the effects 
of gravitational waves in the correlation scalars, however 
we will leave this to future work. Thus, the results pre- 
sented here apply only to scalar perturbations. 

In terms of the scalar perturbations in the cN gauge, 
for a flat FLRW background, the correlation objects (|47|) 
reduce to 



-p(i) = _ 
a 2 



1 



2( «0 2 > + < W-^'f > 



( (VaV b /?') (v a vV) ) 



(54a) 



6( d A ^d A 4> ) + {d A { v - 4>)d A { v - V) } 



(V A \7 B VcP)(V A V B V c f3) ) 



, (54b) 



7?U) +p( 2 ) = 1 



-2H{{v'ip)-{ip'ip) + {ij'{<p-<ip)) 

+( W - V') } + ( (V A V B (3)(V A W B (3') } } 

(54c) 



( d A fi" (8 AV - Hd A f3') ) 



(54d) 
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where f3 is defined in Eqn. (|51b|) . 

Since we are working with a flat FLRW background, 
it becomes convenient to transform our expressions in 
terms of Fourier space variables. This will also high- 
light the problem with large scale fluctuations which was 
mentioned in Sec. IIIII We will use the following Fourier 
transform conventions : For any scalar function f(rj,x), 
its Fourier transform satisfies 

/(**) = / ^ e **/sfo). 

f s {r,)= JcPxe-^ffaO). (55) 

Consider an average of a generic quadratic product of two 
scalars f^(x) and f^ 2 \x) where we have suppressed the 
time dependence since it simply goes along for a ride. 
Using the definition (|4ip . and keeping in mind that the 
scalars are real, it is easy to show that we have 

( / (1 ¥ 2 > >(*) = / ^^wak - M/^f , 

(56) 

where Wi,(fe,af) is the Fourier transform of the window 
function Wl {x, y) on the variable y , and the asterisk de- 
notes a complex conjugate. 

In the present context, the functions / ^ and will 
typically be derived in terms of the initial random fluctu- 
ations in the metric <pt- which are assumed to be drawn 
from a statistically homogeneous and isotropic Gaussian 
distribution with some given power spectrum. In order 
to ultimately make contact with observations, it seems 
necessary to perform a formal ensemble average over all 
possible realisations of this initial distribution of fluctu- 
ations. The statistical homogeneity and isotropy of the 
initial distribution implies that the functions /W and 
f( 2 ) will satisfy a relation of the type 

[ /^C* U = (^) 3 t {3) (h - k 2 )P flh (\ki\) , (57) 

for some function Pf 1 f 2 (k,r]) which is derivable in terms 
of the initial power spectrum of metric fluctuations, and 
where [ ... ] ens denotes an ensemble average and 5^(k) 
is the Dirac delta distribution. 

Applying an ensemble average to Eqn. (|56|) introduces 
a Dirac delta which forces k\ = k%. Further, the normal- 
isation condition on the window function in Eqn. (1411) 
implies that we have 

W L (k = 0,£) = l, (58) 

which means that all dependence on the averaging scale 
and domain drops out, and we are left with 

[ ( / (1) / (2) ) U = / (0^( fc )- (59) 

Note however, that the right hand side of Eqn. (|59p is 
precisely what we would have obtained, had we treated 



the spatial average ( ... } to be the ensemble average 
[ ... ] ens to begin with. Therefore for all practical pur- 
poses, we are justified in replacing all the spatial averages 
in the expressions for the correlation scalars (f54")l . by en- 
semble averages. 

It is convenient to define the transfer function ^k(r]) 
via the relation 

?jgfa) ( 6 °) 

For the calculations in this paper, we shall assume that 
the cN gauge scalars ip(r),x*) and ip{rj,x) are equal 

<p(r), x ) = ip(r), x) , (61) 

a choice which is valid in first order PT when anisotropic 
stresses are negligible (see Ref. [13]). This simplifies 
many of the expressions we are dealing with. The Fourier 
transform of /3 can be written, using Eqns. (|51b|) and 
([BT]). as 

/W = J^gfa)- (62) 

Finally, in terms of the transfer function $/ c ('7) and the 
initial power spectrum of metric fluctuations defined by 

[ miV>h = ( 2 3 * (3) (£ " k 2 )P vl {ki) , (63) 

the correlation scalars (I54[) can be written as (compare 
Eqns (58)-(61) of Ref. (l(|) 



! / ^P v% {k)^ k )\ (64a) 
S^ = -lJ^k 2 P vi (k)(k^l), (64b) 

P w + v m =-^J £s#Pn(k) (* k *' k ) , (64c) 

* (2) = ~ 2 J £#p«wn (** - f *t) • (64d) 

These expressions highlight the problem of having a fi- 
nite amplitude for fluctuations at arbitrarily large length 
scales (k — > 0), which was mentioned in Sec. IIIII As 
a concrete example, consider the frequently discussed 
Harrison-Zel'dovich scale invariant spectrum [34| which 
satisfies the condition 

k 3 P V i(k) = constant . (65) 

Eqns. (|64"1) now show that if the transfer function $fc(rj) 
has a finite time derivative at large scales (as it does in 
the standard scenarios - see the next section), then the 
correlation objects and all diverge due to 
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contributions from the k — > regime. This demonstrates 
the importance of having an initial power spectrum in 
which the amplitude dies down sufficiently rapidly on 
large length scales (which is a known issue, see Ref. [33lj ) . 
Perturbation theory cannot adequately describe the be- 
haviour of inhomogeneities with arbitrarily large length 
scales 35]. Keeping this in mind, we shall concentrate 
on init ial p ower spectra which display a long wavelength 
cutoff [361. Models of inflation leading to such power 
spectra have been discussed in the literature (37[, and 
more encouragingly, analyses of WMAP data seem to in- 
dicate that such a cutoff in the initial power spectrum is 
in fact realised in the universe [381 ] . 

A final comment before proceeding to explicit calcu- 
lations : In addition to picking up nontrivial correlation 
corrections in the cosmological equations, the averaging 
formalism also requires that the "cross-correlation" con- 
straints in Eqns. (fT6|) be satisfied. In the absence of 
vector and tensor modes, it is straightforward to show 
that the statistical homogeneity and isotropy of the scalar 
metric fluctuations implies that these constraints are 
identically satisfied. It will be an interesting excercise 
to analyse the conditions imposed by these constraints 
in the presence of tensor modes; we leave this for future 
work. 



WORKED OUT EXAMPLES 



nos. At the "zeroth iteration" (see Sec. Ill CI) we have 



i£) -*<.)- as 



(66) 



where Ho is the standard Hubble constant, the scale fac- 
tor is normalised so that o(to) = 1, and TL and H are 
related by 



H(a) = aH{a) 



(67) 



The comoving wavenumber corresponding to the scale 
which enters at the matter radiation equality epoch, is 
given by 



k eq — a eq H(a eq ) = H 



20 2 



1/2 



H ■ 10 5/2 , (68) 



where we have set (see Refs. [30, |39| for details) 



photon 



photon 



4.15 x 10 _b /i 



-51,-2 



8 V 11 



(69) 



where h is the dimensionless Hubble parameter defined 
by Hq = lOO/i km/s/Mpc. For all calculations we shall 
set h = 0.72 go]. 



We will now turn to some explicit calculations of the cor- 
relation integrals (backreaction) , which will show that 
the magnitude of the effect remains negligibly small for 
most of the evolution duration in which linear PT is valid. 
At early times, linear PT is valid at practically all scales 
including the smallest scales at which we wish to ap- 
ply General Relativity. As matter fluctuations grow, the 
small length scales progressively approach nonlinearity, 
and linear PT breaks down at these scales. As we will see, 
however, by the time a particular length scale becomes 
nonlinear, its contribution to the amplitude of the met- 
ric fluctuations correspondingly becomes negligible. In 
practice therefore, one can extend the linear calculation 
well into the matter dominated era, with the expectation 
that the order of magnitude of the various integrals will 
not change significantly due to nonlinear effects (see also 
the discussion in the last section). 

The model we will use is the standard Cold Dark Mat- 
ter (sCDM) model consisting of radiation and CDM [30] • 
We will neglect the contribution of baryons, and at the 
end we shall discuss the effects this may have on the final 
results. We shall also discuss, without explicit calcula- 
tion, the effects which the introduction of a cosmological 
constant is likely to have. In the following, f2 r and Q m 
denote the density parameters of radiation and CDM re- 
spectively at the present epoch tq, with r denoting cosmic 
time. fi r is assumed to contain contributions from pho- 
tons and 3 species of massless, out-of-equilibrium neutri- 



A. EdS background and non-evolving potentials 

Before dealing with the full model (which requires a 
numerical evolution) let us consider the simpler situ- 
ation, described by an Einstein-deSitter (EdS) back- 
ground, with negligible radiation and a nonevolving po- 
tential <p = tp(x) (which is a consistent solution of the 
Einstein equations in the sCDM model at least at sub- 
horizon scales at late times [lO]). Although not fully 
accurate, this example requires some very simple inte- 
grals and will help to give us a feel for the structure and 
magnitude of the backreaction. 

With a constant potential, the only correlation object 
which survives is which evolves like ~ a~ 2 , where 

the scale factor refers to the "zeroth iteration" . The con- 
stant of proportionality can be written in terms of the 
BBKS transfer function TBBKs{k/k eq ) [3(J EH, to give 



s{1> - 4 / s^W. 



BBKS 



(k/k, 



(70) 



where we have 1411 



T, 



BBKS 



(x) 



In [1 + 0.171a;] 
(0.171a;) 



1 + 0.284a; 



+ (1.18a;) 2 + (0.399a:) 3 + (0.490a;) 4 



-0.25 



(71) 
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where x = (k/k eq ). 

The integral in Eqn. (|70|) is well-behaved even in the 
presence of power at arbitrarily large scales, for a (nearly) 
scale invariant spectrum. Since we are only looking for an 
estimate, we shall evaluate the integral in the absence of 
a large scale cutoff, and leave a more accurate calculation 
for the next subsection. For the initial spectrum given 
by 



k P<pi {k) 
2ir 2 



A(k/HoY 



(72) 



where the scalar spectral index n s is close to unity, the 
integral in Eqn. (|70|) can be easily performed numerically 
and has the order of magnitude 



dk 
2^2 



k A P V i(k)Tg BKS 



(k/k eq ) r*j A (k e q) 



AH* 



10 5 , 
(73) 

upto a numerical prefactor of order 1. Since the ampli- 
tude of the power spectrum is A ~ 10 -9 [42j], the overall 
contribution of the backreaction is 



H 2 



1 



(74) 



Now, as long as the correlation objects give a negligible 
backreaction to the usual background quantities, when 
we proceed with the next iteration, the effect of the back- 
reaction on the evolution of the perturbations will also 
remain negligible (at least at the leading order). Hence 
in practice there will be essentially no difference between 
the zeroth iteration and first iteration perturbation func- 
tions. This amounts to saying that when the backre- 
action is negligible, convergence to the "true" solution 
for the scale factor at the leading order, is essentially 
achieved in a single calculation. 



B. Radiation and CDM without baryons 

Let us now turn to the full sCDM model (without 
baryons). An analytical discussion of this model in vari- 
ous regions of (k, ?y)-space, can be found e.g. in Ref. [3Cj |. 
Since we are interested in integrals over k across a range 
of epochs n, it is most convenient to solve this model nu- 
merically. It is further convenient to use (In a) in place 
of r], as the variable with which to advance the solution. 
Also, it is useful to introduce transfer functions like 
for all the relevant perturbation functions in exactly the 
same manner (see Eqn. (1601) ). namely by pulling out a 
factor of ip-j--, since the initial conditions are completely 
specified by the initial metric perturbation. For a generic 
perturbation function s^(-q) (other than the metric fluc- 
tuation ifj:) the transfer function corresponding to s will 
be denoted by a caret, so that 



(75) 



The relevant Einstein equations can be brought to the fol- 
lowing closed set of first order ordinary differential equa- 
tions (adapted from Eqns. (7.11)-(7.15) of Ref. [20]), 



ggfc 
d(lna) 



dS k 
d(lna) 



K 2 \ 
3E 2 J 



In 



I ( A 4 . 
-——j— fl m Sk H — r 6os; 
2E z a \ a 



K, 
E 



V k +3 



d(lna) 



^e 2L = _^e , 

dQna) E xk d(bxa) 



K 



d(lna) 3E 



Ok 



dV k 
d(lna) 



-Vk + 



(76a) 



(76b) 



(76c) 



(76d) 



(76e) 



Here we have introduced the dimensionless variables 



K^— ; E(a) 



H(a) H(a) 



Ho 



Ho 



(77) 



and the various perturbation functions are defined as fol- 
lows : 8k is the fc-space density contrast of CDM, 6ofc and 
Qik are the monopole and dipole moments respectively 
of the fe-space temperature fluctuation of radiation, and 
(— iVk) is the fc-space peculiar velocity scalar potential of 
CDM (i.e., the real space peculiar velocity is va = 9av 
where v is the Fourier transform of (— iVk))- 

Assuming adiabatic perturbations, the initial condi- 
tions satisfied by the transfer functions at a — ai are 
(adapted from Ch. 6 of Ref. [U) 



$k(o») = 1 ; S k (ai) 



V k (ai) = 38ifc(ai) = - 



3 
2 

1 K 



<dok(ai 



1 



2 E(ai) 



(78) 



We choose = 10 -16 , which corresponds to an ini- 
tial background radiation temperature of T ~ 10 3 GeV. 
While this is not as far back in the past as the energy scale 
of inflation (which is closer to ~ 10 15 Gev), it is on the 
edge of the energy scale where known physics begins [33| . 
This makes Eqn. (|69p unrealistic since we have ignored 
all of Big Bang Nucleosynthesis and also the fact that 
neutrinos were in equilibrium with other species at tem- 
peratures higher than about lMev. However the modifi- 
cations due to these additional details are not expected 
to drastically change the final results, and these assump- 
tions lead to some simplifications in the code used. The 
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FIG. 1: The evolution of two large scale modes. Also shown 
is the Kodama-Sasaki analytical solution in the large scale 
limit kn < 1, given by Eqn. ([79]) . 



FIG. 2: The transfer function normalised by its constant 
value at large scales, at the epoch a — 500a e9 . The dotted 
line is the BBKS transfer function (171|l . 



goal here is only to demonstrate an application of the for- 
malism; more realistic calculations accounting for the ef- 
fects of baryons can also be performed (see, e.g. Behrend 
et al. [l(J who incorporate these effects for the post- 
recombination era, albeit in the Buchert formalism). 

In order to partially account for the fact that infla- 
tionary initial conditions are actually set much earlier 
than a = 10 -16 , we impose an absolute small wavelength 
cutoff at the scale which enters the horizon at the ini- 
tial epoch which we have chosen. In the above notation 
this corresponds to setting K max = E(ai) ~ 10 13 . This 
makes sense since scales satisfying K ^> K ma x have al- 
ready entered the horizon and decayed considerably by 
the epoch a = 10~ 16 . There is a source of error due to 
ignoring scales K > K max which have not yet decayed 
significantly, but this error rapidly decreases with time as 
progressively larger length scales enter the horizon and 
decay. [In fact, in practice to compute the integrals at 
any given epoch a — a*, one only needs to have followed 
the evolution of modes with K < ~ 5000E(a*) : more 
on this in the next subsection.] More important is the 
cutoff at long wavelengths, which we set at K m i n = 1 
(corresponding to fc m j n = Hq), which is firstly a natu- 
ral choice given that Hq 1 is the only large scale in the 
system, and is secondly also guided by analyses of CMB 
data which have detected such a cutoff [H| . We will see 
that reducing K m i n even by a few orders of magnitude, 
does not affect the final qualitative results significantly. 



1. Numerical Results 

Equations ([75]) with initial conditions ([75)) were solved 
with a standard 4th order Runge-Kutta integrator with 
adaptive stepsize control (based on the algorithm given 
in Ref. [IH). For the integrals in Eqns. ([M]) . only the 
function < &fc(a) needs to be tracked accurately. Hence, 



although Gofe and Oi^ are difficult to follow accurately 
beyond the matter radiation equality a eq = (Q r /f2 m ) ~ 
8 x 10~ 5 due to rapid oscillations, the integrals can still be 
reliably computed since Qok and Gtu- do not significantly 
affect the evolution of $^ in the matter dominated era 
(as seen in Eqn. (|76a|) ). 

To see that known results are being reproduced by the 
code, consider Figs. [1] and [21 as examples. Fig. Q] shows 
the evolution of two scales corresponding to K = 1 (k = 
Ho Mpc^ 1 ) and K = 0.01. The first enters the horizon at 
the present epoch, while the second remains superhorizon 
for the entire evolution, satisfying kr) <C 1. In this limit 
an analytical solution exists in the sCDM model, due to 
Kodama and Sasaki [3(| EH , given by 



*fc(») 



1 



10y 3 



16-v/l + y + 9y 3 + 2y 2 - 8y - 16 



where y = a/a eq , and this function is also shown. 



(79) 
Clearly 

all the curves in Fig. [T]are practically identical. 

Fig. [2] shows the function $fc normalised by its (con- 
stant) value at large scales, at the epoch a = 500a eg ~ 
0.04, (which is well into the matter dominated era). The 
dotted line is the BBKS fitting form given in Eqn. (|71p 
with k eq given by Eqn. (|68j) . 

To numerically estimate the integrals in Eqns. ([64| . 
the values of <E>fc and its first and second derivatives with 
respect to (In a) are needed across a range of K values. 
For reference, note that the following relations hold for a 
generic function of time w(r]), 



drj da d(\n a) 



(80) 



Based on the earlier discussion, the initial power spec- 
trum P(pi(k) is taken to satisfy 



k^ P(pi (&) 



A , for H Q < k < km 



W(ai), (81) 



14 



le+09 
le+08 - 
le+07 
le+06 
le+05 - 
le+04 
1000 
100 - 
10 
1 

0.1 - 



- "eq 

= 0.0125a„ 



0.01 

0.0001 0.001 




0.01 0.1 1 10 100 1000 10000 

k(h Mpc" 1 ) 



0.1 



0.001 



"eq 

a = 0.0125a 



le-05 

0.0001 0.001 0.01 0.1 1 10 100 1000 10000 

k (h Mpc" 1 ) 



FIG. 3: The dimensionless integrand of S (1 \ namely the 
function (k/ Hq) 2 <&\, at three sample values of the scale 
factor. The function dies down rapidly for large k, with the 
value at some k being progressively smaller with increasing 
scale factor. The declining behaviour of the curves for 
a = a eq and a = 200a e<3 extrapolates to large k. 



and zero otherwise, and we set 

A = 1.0 x 1(T 9 . 



(82) 



which, for the sCDM model follows from the convention 
(see Eqn.(6.100) of Ref. Op) A = (5fe/3) 2 with S H ~ 
2 x 1(T 5 (see, e.g. Ref. (Hf). 

Consider Figs. [3]and[2J which highlight two issues dis- 
cussed earlier. Fig. [3]shows the integrand of S^ 1 ' at three 
sample epochs, and we see that the integrand dies down 
rapidly at increasingly smaller k values for progressively 
later epochs. (The other integrands, not displayed here, 
also show this rapid decline for large k.) [We have not 
shown the integrand at the later two epochs for all val- 
ues of k since this was computationally expensive, but 
the declining trend of the curves can be extrapolated to 
large k, which is well understood analytically [3fJ.] This 
justifies the statement in the beginning of this section, 
that at any epoch a* it is sufficient to have followed the 
evolution of scales satisfying K < 5000i?(a*) for comput- 
ing the integrals. Secondly, Fig. [4] shows the behaviour 
of k z / 2 \5k\ = A x / 2 \5k\ at the same three epochs, and com- 
paring with Fig. [3]we see that at any epoch, the region of 
/c-space where linear PT has broken down, does not con- 
tribute significantly to the integrals. This is in line with 
the conjecture in Ref. [45[ that the effects of the backre- 
action should remain small since the mass contained in 
the nonlinear scales is subdominant. 

Due to the structure of the integrals and the chosen 
initial power spectrum, it is convenient to compute the 
integrands in Eqns. (|64p equally spaced in (InK), and 
then perform the integrals using the extended Simpson's 



rule [43j. If 2 N + 1 points are used to evaluate a given 
integral, resulting in a value In say, then the error can 
be estimated by computing the integral with 2 N ~ 1 + 1 
points to get Ziv-i, and estimating the relative error as 



FIG. 4: The dimensionless CDM density contrast. Together 
with Fig. [3] this shows that nonlinear scales do not impact 
the backreaction integrals significantly. 



\Tn-i/Zn\ — 1- With N = 10, the estimated errors in all 
the integrals at all epochs were typically less than 0.1%. 
A bigger error is incurred in computing the integrand 
itself at any given epoch, leading to estimated errors of 
order ~ 1% in V (1) and T {1) + V (2 \ with a larger 
error in iS*- 2 ) as explained below. 

The second derivative d 2 $k/d(lna) 2 proves to be dif- 
ficult to track numerically. At early times, when most 
scales are superhorizon, the Kodama-Sasaki analytical 
solution ([75)) is a good approximation for most values of 
k. Using this one can see that at early times the value 
of the derivative is numerically very small, and is dif- 
ficult to reliably estimate due to roundoff errors. For 
this reason the integral <S^ 2 ' could not be accurately es- 
timated at early times. However, the structure of the 
integrand of (|64d[) shows that the largest contribu- 
tion comes from large (superhorizon) scales (the small 
scales being subdominant due to the presence of $fc 
and 1/fc 2 ). An analysis using the Kodama-Sasaki solu- 
tion then shows in a fairly straightforward manner that 
the behaviour of the backreaction term is /H 2 \ ~ 
10~ 6 (a/a eq )(Ho/k m i n ) 2 for our choices of parameters, 
where (a/a eq ) <C 1. At intermediate times around 



l eq 



and later, although it becomes computation- 
ally expensive to obtain conver gen t values for the second 
derivative at all relevant scales [46| , moderately good ac- 
curacy (1-5%) can be achieved. 

The results are shown in Fig. [51 in which the magni- 
tudes of the correlation integrals of Eqn. ([tj4"[) . normalised 
by the Hubble parameter squared H 2 (a) = (TL/a) 2 are 
plotted as a function of the scale factor in a log- log plot. 
The values for are shown only for epochs later than 
a ~ 0.01a e? ~ 10~ 6 . We see that at all epochs, the 
correlation terms remain negligible compared to the cho- 
sen zeroth iteration background. Also, in the radiation 
dominated epoch all the correlation scalars (except 
whose evolution couldn't be accurately obtained) track 
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FIG. 5: The correlation scalars ( "backreaction" ) for the 
sCDM model, normalised by H 2 (a). S w , V (1) and <S (2) are 
negative definite and their magnitudes have been plotted. 
The vertical line marks the epoch of matter radiation 
equality a = a eq . 

the ~ a -4 behaviour of the background radiation den- 
sity (see also the second paper in Ref. [|J). The dis- 
cussion above shows however that the magnitude of 
is far smaller than the other backreaction functions at 
early times, for a cutoff at k m i n = Hq. On the other 
hand, in the matter dominated epoch iS' 1 ' dominates the 
backreaction and settles into a curvature-like ~ a~ 2 be- 
haviour (note that in the matter dominated epoch we 
have H 2 ~ a~ 3 ). This can also be compared with the 
results of Ref. [3j. As for the signs of the correlations, 
«S« S< 2 >, and are negative throughout the evolu- 
tion while + V^ 1 is positive throughout. 

Finally, a few comments regarding the effects of ig- 
noring baryons, nonlinear corrections, etc. Including 
baryons in the problem (with a background density pa- 
rameter of fib ~ 0.05) will lead to a significant suppres- 
sion of small scale power (by introducing pressure terms 
which will tend to wipe out inhomogeneities) and also a 
small suppression of large scale power. This effect causes 
a (downward) chan ge i n the late time transfer function 
of roughly 15-20% [30], and therefore cannot increase 
the contribution of the backreaction. Quasi-linear cor- 
rections can lead to significant changes in the transfer 
function, but do not cause shifts by several orders of 
magnitude (see Ref. [13] and references therein). Hence 
accounting for changes due to quasi-linear behaviour will 
also not increase the magnitude of the backreaction by 
a large amount (see also Ref. [13]). As for effects from 
fully nonlinear scales, we have seen that these can be 
expected to remain small, or at least not orders of mag- 
nitude larger than those from linear scales (see also the 
discussion in the last section, and Ref. [H). 

Adding a cosmological constant (and retaining a flat 
background geometry) will change the qualitative fea- 
tures of the correlation functions by shifting the scale 
k eq (due to a reduced fi m , which will also increase the 



power spectrum amplitude [3fjj . but again not by orders 
of magnitude). Also, the late time behaviour of the cor- 
relation scalars will be affected since the potential $fc 
will decay at late times instead of remaining constant. 
Regardless, the backreaction is expected to remain small 
even in this case (which is also indicated by the calcu- 
lations of Behrend et al. [HI in the Buchert framework 
0). 



VI. DISCUSSION 

This paper has presented an analysis of cosmological 
perturbation theory (PT) in the fully covariant aver- 
aging framework of Zalaletdinov's Macroscopic Gr avity 
(MG) [1, 0] and its restriction to spatial averaging [2lf . 
While this framework is generally covariant, the issue of 
gauge dependence in perturbation theory introduces cer- 
tain subtleties in the problem. Wc have shown that, pro- 
vided one takes seriously the idea that the cosmological 
background must be defined by an averaging procedure 
[23| . it is possible to attach a gauge invariant meaning 
to the averaging condition and the corresponding cor- 
relation objects which appear as corrections to the cos- 
mological equations. While there remains considerable 
freedom in an explicit choice of the averaging operator 
(through a choice of the volume preserving gauge used in 
its definition), this freedom can be fixed by some addi- 
tional requirements based on knowledge of cosmological 
PT in the standard implementations. In particular we 
have seen that properties of the conformal Newtonian or 
Poisson gauge can be used to motivate a fully specified 
choice of the averaging operation adapted to first order 
PT. 

One prerequisite to the formulation of a consistent av- 
eraging framework in the context of perturbation theory, 
is the absence of perturbative fluctuations with arbitrar- 
ily large wavelengths, since such fluctuations would ren- 
der meaningless the notion of recovering a homogeneous 
background on averaging. This problem also manifested 
itself in the correlation integrals (|64p . which diverge in 
the presence of a finite amplitude of fluctuations as the 
wavenumber k — ► 0. Accordingly, all the calculations 
of this paper have assumed that the initial power spec- 
trum of metric fluctuations has a sharp cutoff at the scale 
corresponding to k — a hypothesis which is in fact 
supported by analysis of CMB data [38j . 

The main purpose of this paper was to lay down the 
formalism of MG in a language most convenient from the 
point of view of cosmological PT. This was accomplished 
by writing Eqns. (|47|) , (|54|) and the Fourier space ver- 
sion (|64[) (with certain simplifying assumptions regarding 
vector and tensor perturbations which can if needed be 
relaxed in a completely straightforward manner). This 
was supplemented by calculations in the sCDM model 
[3(i| (which is the flat FLRW model with radiation and 
Cold Dark Matter but no Dark Energy) with the addi- 
tional simplification of ignoring the baryons. The ana- 
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lytical results of Sec. IV Al as well as the more detailed 
numerical results of Sec. IV Bl show that the correlation 
objects or "backreaction" remain negligibly small upto 
epochs corresponding to a scale factor of a ~ 0.01. While 
the calculations ignored corrections from quasi-linear and 
nonlinear scales, these are not expected to contribute dra- 
matically to the correlations obtained here, an expecta- 
tion which is justified by the work of Bchrcnd et al. [10] 
and further by the calculation in Ref. [48[ (see below). 

We have seen that by using the framework of MG, we 
have completely bypassed the problem mentioned in the 
Introduction, which one faces when applying the Buchert 
framework to cosmological PT, namely of having to deal 
with two scale factors. Here, one has a single well-defined 
scale factor associated with the background metric, and 
its evolution can be obtained in an iterative fashion as 
described in Sec. Ill CI In practice, we saw that since 
the backreaction is small, convergence can be achieved 
by essentially a single calculation, at least in the context 
of first order PT. 

This brings us to a final, and very important issue : 
What is the magnitude and behaviour of the backreac- 
tion in the fully nonlinear regime of structure formation? 
There are (at least) two possible avenues to approach this 
question. The first is to set up the problem in a manner 
which is suitable for TV-body simulations. The iterative 
approach suggested earlier can presumably be adapted 
to full-fledged iV-body codes as well. While this is a 
possibility worth pursuing, there is also a less involved 
(but correspondingly less realistic) way of determining 
the effect of nonlinearities on the backreaction, which is 
to study toy models of structure formation. As men- 
tioned earlier, such a toy model of spherical collapse was 
recently presented in Ref. [20| . and it was shown by an 
explicit coordinate transformation, that the metric can 
be brought to the Newtonian form Jl}. Now, it is worth 



noting that while all the calculations of this paper as- 
sumed that first order PT is valid, the actual expressions 
for the correlation integrals in Eqns. (|54p only assume 
that the potentials ip and yb satisfy |<p|,|V'l "C 1. The 
dynamics governing the potentials is irrelevant at this 
stage. This means that, as long as one is interested in 
leading order effects only, the expressions in Eqns. (|54|) 
can be directly applied to any model of structure forma- 
tion where the metric can be brought to the conformal 
Newtonian form [49|: in particular they can be applied to 
the model of Refijll- This has been done in Ref. [II], 
and one finds that even in the fully nonlinear regime, 
the effect of the backreaction remains negligible. In this 
context see also Ref. 0] for a third approach. 

To conclude, all our calculations and arguments seem 
to indicate that the averaging of perturbative inhomo- 
geneities in a consistent manner appears to lead only to 
very small effects. This does not, however, mean that 
the effects do not exist. It remains to be seen whether 
any observable consequences of the backreaction may be 
detectable by future experiments. 
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